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Modeling of rough frictional interfaces is often based on asperity models, in which the 
behavior of individual microjunctions is assumed. In the absence of local measurements 
at the microjunction scale, quantitative comparison of such models with experiments 
is usually based only on macroscopic quantities, like the total tangential load resisted 
by the interface. Recently however, a new experimental dataset was presented on the 
onset of sliding of rough elastomeric interfaces, which includes local measurements of 
the contact area of the individual microjunctions. Here, we use this more comprehensive 
dataset to test the possibility of quantitatively matching the measurements with a model 
of independent asperities, enriched with experimental information about the area of 
microjunctions and its evolution under shear. We show that, despite using parameter 
values and behavior laws constrained and inspired by experiments, our model does not 
quantitatively match the macroscopic measurements. We discuss the possible origins of 
this failure. 


Keywords: rough contact, elastomer friction, onset of sliding, asperity model, shear-induced area reduction, 
stick-slip, elastic interactions 


1. INTRODUCTION 


The mechanical behavior of contact interfaces between rough solids is crucial to understand their 
tribological properties. The rough contact mechanics community has been developing models in 
two main directions (see Vakis et al., 2018 for a recent review). First, asperity models in which 
the contact interface is divided into well-defined microjunctions actually carrying the normal and 
tangential loads applied to the contacting solids (Braun and Réder, 2002; Ciavarella et al., 2006; 
Violano and Afferrante, 2019). Each microjunction is ascribed a set of individual properties (e.g., 
its height, radius of curvature or friction coefficient) necessary to apply some assumed behavior laws 
[e.g., any contact (Johnson, 1987) or friction law (Le Bot et al., 2019)] when submitted to an external 
stimulus. The macroscopic behavior of the interface is then the emerging, collective response of the 
population of microjunctions (Tromborg et al., 2014; Braun and Peyrard, 2018; Costagliola et al., 
2018). Second, continuum models in which the input quantity is the full topography of the rough 
surfaces, and an exact solution of the unilateral contact and friction problem is seeked (Pastewka 
and Robbins, 2014; Yastrebov et al., 2017; Ponthus et al., 2019), again under some assumptions on 
the interfacial behavior, concerning, e.g., elasticity, friction, and adhesion. 
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Each approach can be used to produce two types of results, in addition to the macroscopic loads on the interface, they 
either deterministic or statistical. Deterministic results are include the evolution under shear of the individual contact areas 
obtained for a given topography (for continuum models) or a and shapes of the many microjunctions forming the interface 
given set of model parameters (for asperity models, including (Figure 1B). 
the properties of each microjunction) and are thus specific to The philosophy of this work is to start with a model of 
those input data. They are relevant for quantitative comparison independent asperities like the earthquake model of Braun 
with a particular experiment. In contrast, statistical results and Peyrard (2008), enrich it with the recently identified 
are the expected results of a large number of deterministic phenomenology of shear-induced area reduction, and genuinely 
calculations performed on statistically similar random surfaces. ask the question whether such a model is sufficient to 
In asperity models, statistical results are obtained when using quantitatively match a particular experimental dataset. In other 
probability density functions (pdfs) of the microjunction words, we do not aim at a definitive model of the incipient 
properties (Greenwood and Williamson, 1966; Braun and tangential loading and onset of sliding of rough elastomer 
Peyrard, 2008; Thogersen et al., 2014). In continuum models, contacts. Rather, we make one single step ahead compared 
they are usually obtained using the power spectrum density (psd) to the models in the literature, and try to conclude whether 
of the topographies under study (Persson, 2001). In the following, this step (including shear-induced area reduction) is sufficient 
we aim at finding a quantitative match with a specific set of | or not. Such an approach can only be fruitful if the values 
measurements, so we will consider deterministic results. of the model parameters are sufficiently constrained by the 

Both asperity and continuum models have been widely experimental dataset, so that one avoids fortuitous agreement. 
explored in the context of rough contacts under purely normal This can be achieved (i) by limiting to the strict minimum 
load, with a recent study explicitly comparing the relative merits the number of parameters that cannot be directly measured 
of the two approaches (Miiser et al., 2017). Several studies experimentally, and (ii) by performing a thorough exploration 
aimed at a quantitative comparison between deterministic model of the parameter space for those remaining, unconstrained 
results and local, microjunction level measurements (see e.g., parameters. In this work, we did our best to apply this strategy, 
McGhee et al., 2017; Acito et al., 2019). In contrast, to our which in our case leads to an unsatisfactory agreement. This 
best knowledge, such comparisons have not been reported in result is nevertheless a progress in the sense that it clarifies the 
the case of sheared multicontacts. Here we will attempt to range of assumptions that remain to be questioned and improved 
build an asperity model able to quantitatively match recent in future studies. 
measurements performed on the incipient tangential loading In section 2.1, we describe the asperity model and provide the 
and onset of sliding of a rough elastomer slab in contact with experimental constraints on the model parameters in section 2.2. 
a smooth glass plate (Sahli et al., 2018, 2019) (Figure 1A). Quantitative comparisons between the model and measurements 
Those measurements (see a typical example in Figures1C,D) are given in section 3, while in section 4 we discuss the possible 
are particularly interesting and constraining for models because, reasons for the absence of good matching between the two. 
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FIGURE 1 | Experiments that our model attempts to reproduce. (A) Sketch of the experimental setup. (B) Typical segmented image of the interface showing individual 
microcontacts in white, for P = 6.40 N. (C) Concurrent time evolutions of the tangential force Q (red) and the area of real contact (blue), for P = 3.10 N. (D) Area of 
real contact as a function of the tangential force, for the same data as in (C). 
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2. MATERIALS AND METHODS 
2.1. Model Description 


We consider the frictional interface between a slider of mass M 
and a track. The tangential displacement of the slider, X(t), is 
assumed to obey the following equation of motion: 





MX + MnX = ky (vt — X) — F, (1) 


where M is the sliders mass, kz, is the stiffness of the 
loading spring through which the slider is pulled at constant 
velocity v, 7 is an effective viscous parameter accounting for 
dissipation, e.g., in the air or in the loading spring, the dot 
indicates the time derivative, and F is the resistive force due to 
interfacial/adhesive friction. 

We assume that the interface is a multicontact made of N 
independent individual microjunctions, each resisting a force fi, 
so that F = ae fi. Each microjunction can be in either of 
two states. First, a pinned state during which the junction acts 
like a (time-evolving) elastic spring of stiffness k;, so that fi = 
ki(X(t) — xi), with x; the slip displacement of the junction with 
respect to the track (e.g., x; = 0 as long as junction i has never 
been slipping). When a threshold force fs; is reached, the junction 
enters a slipping state, during which f; = € fi. Note thate < 1, 
so that fs; and € fs; are the analogs, at the junction level, of a static 
and a dynamic friction force, respectively. 

The mechanical behavior of individual junctions is inspired by 
experimental observations made on the same setup and materials 
in contact as in Figure 1A, but when the rough slab is replaced 
by a single smooth sphere (Sahli et al., 2018, 2019; Mergel 
et al, 2019). The resulting sphere/plane contact is assumed 
to be representative of an individual microjunction within a 
multicontact like that of Figure 1B. Those experiments, carried 
out both for large normal loads (Sahli et al., 2019) and for small 
(even negative) normal loads (Mergel et al., 2019), have shown 
that, under increasing shear, the initially circular contact shrinks 
anisotropically and becomes increasingly ellipse-like. As shown 
in Sahli et al. (2019) the shrinking minor axis of the ellipse is 
parallel to the shear loading direction, while the variations of the 
major axis (in the direction orthogonal to shear) can be neglected. 

Defining £y; and ¢,; the sizes of an elliptic microjunction 
along and orthogonal to shear, respectively, we can define its area 
as Aj = Teil Li. Following Mindlin (1949), the stiffness of such 
an elliptic contact along the shear direction is, assuming no-slip 
contact conditions: 


Te LE 


"= +») [ko — 4(K@ - E(e))] | 


(2) 





with E and v being the Youngs modulus and Poisson’ ratio of 


the material that constitutes the microjunctions, e = ,/1 — +Ë is 


the excentricity of the junction, K and E are the elliptic integrals 
of the first and second type, respectively. Note that assuming that 
microjunctions are elliptic is the simplest increment of realism 
compared to a circular assumption, in order to account for the 
complex shapes observed for microjunctions in the experiments. 


Assuming that each microjunction is initially circular, we can 
define the common initial value, £0;, of £ 1; and ¢\; from its initial 


individual area Ao; as Loi = ,/ Aoi As already mentioned, €]; 


varies negligibly under shear, so we will consider that £1; = Loi 
at all times. The evolution of each £j; is then deduced from the 
shear-induced area reduction reported in Sahli et al. (2018): 


1 
Aj = Aoi — ap? (3) 
Ai 


with œp and the exponent p two constant parameters of the 
model. The size of junction i along the shear direction is 
thus simply ¢); = Hi, Note that there is currently no 
rigourous contact mechanics theory for the evolution of the shear 
stiffness of a sheared sphere/plane contact that would incorporate 
anisotropic contact area reduction. Here, such a behavior is 
approximated at all times by the combination of Equation (2), 
which is valid under no-slip assumption, and of Equation (3), 
which was empirically found at macroscale. Doing so, we assume 
that Equation (3) also applies at microscale, as suggested by the 
existence of common values of œ, and p for both the macro- and 
micro-scales (Sahli et al., 2018). 

For each microjunction, Equation (3) is used from the 
beginning of the experiment, when f; assumed to be 0, up to when 
the junction first starts to slip (when f; = fsi). At that instant, A; 
takes the value As; = Aoi — @p T J- For later times, based on the 

0i 


observation of the typical behavior of A; during the experiments 
of Sahli et al. (2018) (see Figure 2), we assume that A; always 
remains equal to Asi. 

In contrast, the force resisted by a microjunction can vary with 
time after the first onset of slipping. When the slider’s velocity, 
X(t), gets smaller than a minimum value Xmin = Cmin X v, With 
Cmin à scalar parameter, we assume that all the slipping contacts 
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FIGURE 2 | Left axis: time evolution of the contact area A; of five typical 
microjunctions in the experiment at P = 4.01 N. Right axis: concurrent 
evolution of the macroscopic tangential force Q. 
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will repin, with a position x; = X—e€fsi/ki. Doing so, at repinning, 
there is no force discontinuity, as the repinning force k;(X — xi) 
is equal to the slipping force efsi. 

Following Sahli et al. (2018), the threshold force fs; at which 
the microjunction starts to slip is assumed to be proportional to 
its area at the same instant, i.e., fsi = o Asi, with o the frictional 
shear strength of the contact. 

The algorithm used to solve the model numerically is provided 
in Appendix 1. 


2.2. Experimental Constraints 

In order to quantitatively compare the model results with the 
multicontact experiments reported in Sahli et al. (2018, 2019), 
we need to feed the model with parameter values based on the 
measurements. In the following list, we first provide all constant 
parameter values that are directly accessible experimentally, with 
the error bars when relevant. 


e M = Mo + Mi, where Mo = 100 g is the mass of the slider and 
Mi = 0,55, 111, 215, 308, 552 g an additional mass, for the six 
experiments performed. All masses are given at +1g. 

e kr = 9,200 + 200 N/m (Sahliet al., 2018). 

e v=0.1 mm/s. 

e E = 1.6 + 0.1 MPa (Sahli et al., 2018). This value and the 
associated error bar are the mean value and standard deviation 
over 32 estimates using 5 different spherical PDMS samples 
prepared in the same conditions. 

e v is assumed to be equal to 0.5, as is classically done 
for elastomers. 

e the individual values of the initial areas of all microjunctions, 
Aoi, are extracted from the initial image (for Q = 0), 
segmented as described in Sahli et al. (2018). The fact that 
all microjunctions have different areas is the result of the 
random nature of the elastomer surface topography (see a 
typical Power Spectrum Density in Supplemental Material of 
reference Sahli et al., 2019) and of its elastic contact interaction 
with the rigid glass plate. 

e N is also extracted form the same segmented image. 

© Oexp = 0.23 + 0.02 MPa (Sahli et al., 2018), is the experimental 
value of the frictional shear strength of the interface, 
determined from a linear fit of (As, Qs) for the 6 experiments. 
Q; is the macroscopic static friction (peak) force and À, is the 
total area of real contact at the same instant. We will discuss 
below how the value of ø in the model is related to ep. 








There are three model parameters which cannot be directly 
measured in experiments: 7, €, and Cmin. 

n is introduced to enable energy dissipation in the system, thus 
avoiding spurious oscillations of the slider. However, the value of 
n should not be too large, because it would prevent the possibility 
of stick-slip in the model. We found that stick-slip exists up to 
n between 180 and 200, but for those large values, the initial 
stick-slip cycles are significantly different from the experiments. 
In practice, we found that 


n = 100 (4) 


is a good compromise between oscillation reduction and a 
reasonable reproduction of the stick-slip sequence. The results 
are rather insensitive to the precise value of n, since n = 50 gives 
virtually identical results. 

€ has a leading order control on the amplitude and period of 
the tangential force fluctuations during stick-slip. Systematic tests 
of the model for various values of € led us to choose 


€ = 0.90. (5) 


In particular, this value is sufficiently small to enable stick- 
slip for all six normal loads (as observed experimentally), while 
reproducing reasonably well the amplitude and period of the 
stick-slip sequences in all cases. 

Note that in the model, if there was no stick-slip, the 
steady-state sliding friction force would be pad 1€0Asi (all 
microjunctions are in their slipping state). In order for this value 
to match the macroscopically measured value Q; = OexpAs, one 
has to impose that 

ome (6) 
€ 
and this is what we do in the following. 

Our tests showed that the value of cm, has no impact on the 
results as far as it is sufficiently small. For instance, simulations 
with Cmin =10~> are essentially undistinguishable from those 
using 0.01. The reason is that, when |X(t)| crosses the value 
Cmin X v, the velocity drop is so fast that the time at which the 
crossing occurs is almost independent on the value of cyin. In 
our calculations, we will use 


Cmin = 0.01. (7) 


Extracting values for p and a, in Equation (3) requires fitting 
the power law relationship between the individual area reduction 
parameters, œ;, and the initial areas, Ao;, presented as purple 
squares in Figure 3 of Sahli et al. (2018). Such a fitting is actually 
difficult due to the large dispersion of the data, as can be inferred 
from the large difference in total area decay of microjunctions 1- 
3 in Figure 2, although they start with almost identical areas. A fit 
letting both œ, and p as fitting parameters gives 95% confidence 
error bars as large as 600% for the optimum value for œp, which 
is not a viable option. We then tried to fix the value of p and fit 
the data with a, being the only fitting parameter. We found that 
the quality of the fit (quantified by its R? value) was essentially 
independent of p (as long as it is not too different from the 
value 3/2 proposed in Sahli et al., 2018), preventing any objective 
choice of p. 

Based on those observations regarding the determination of p 
and a, from experimental data, in our model studies we decided 
to fix p and, for each value of p, we determined the value of œp 
that gives the best agreement between the area decay predicted 
by the model and that measured in the experiments. To do that, 
we fitted both the experimental and model version of the curve 
A(Q) by a quadratic function of the form A = Ag — aQ. 
Ao being the same in the model as in the experiment (because 
Ao = - Aoi), the fitting procedure enables identification of 
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an œp which provides an exact match between the two quadratic 
decays. Importantly, we found that, for all tested values of p 
close to 3/2 (the value suggested in Sahli et al., 2018), the model 
results (when using the corresponding fitted œp) were almost 
undistinguishable. So, in practice, we chose p = 3/2, for which 
the model studies give an optimal a, = 0.45 1071 m°/N? for 
the experiment with the smallest normal load, and a, = 1.00 
101$ m°/N? for the experiment with the largest normal load. We 
then adopted the average value between both, œp = 0.725 10715 
m°/N?, as a constant to be used for all experiments. 


3. RESULTS: QUANTITATIVE COMPARISON 


We run the model of section 2.1 with the parameter values 
described in section 2.2, for the six different PDMS/glass 
multicontact experiments reported in Sahli et al. (2018). Figure 3 
compares, for two different normal loads, the measured time 
evolution of the area of real contact and tangential force to 
their corresponding model predictions. Note that the initial real 
contact area is essentially proportional to the normal load, as 
widely discussed in the contact mechanics modeling literature 
[see e.g., the reviews (Persson et al., 2005; Vakis et al., 2018)] 
and confirmed in the experiments discussed here (see Figure S2 
of Sahli et al., 2018). To facilitate comparison between model 
prediction and measurement, the time origin of the experimental 
data has been offset by the amount necessary to superimpose the 
measured and predicted force curves in the central portion of 
their initial increase. Note that the initial non-linear increase of 
the measured force is due to the non-vanishing bending stiffness 
of the steel wire used to pull the slider, when it first bends around 
a pulley before a significant tension arises along the wire. The 
apparent difference between the measured and predicted values 
of the initial area of real contact is due to the above mentioned 
time offset: the initial predicted value exactly corresponds to 
the measured value from the first image, but the latter image 
now corresponds to a negative time and is thus not shown in 
the figure. The observed difference is of the order of the area 
measurement noise, presumably due to temporal fluctuations in 
the illumination and noise in the cameras sensor. 

Figure 4 then shows, for all normal loads, the evolution of 
the area of real contact as a function of the tangential force, for 
both the measured and predicted data. This figure is similar to 
Figure 2A in Sahli et al. (2018), but shows all measurements 
points rather than just 1 of 130. Note that stick-slip is responsible 
for the accumulation of nearly horizontal cycles close to the 
minimum area/maximum force point of each curve. Also note 
that the model forces can transiently exceed the value 0:,,A, but 


exp 
= A, as expected. 





always remain smaller than o A = 


4. DISCUSSION 


Although other combinations of model ingredients may have 
been proposed, we believe that our model incorporates all of 
the currently available knowledge on the system that we tried to 
reproduce. As such, it can be seen as the most comprehensive 
independent asperities model of shear multicontacts so far, to be 
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FIGURE 3 | Direct comparison between measurements and model 
predictions, for two typical experiments with either P = 1.53 N (A) or P = 3.10 
N (B). Time evolutions of the measured (dashed, black) and predicted (solid, 
blue) area of real contact, and of the measured (dashed, magenta) and 
predicted (solid, red) tangential force. 











used for deterministic comparison with the experiments of Sahli 
et al. (2018, 2019). 

Most of the model parameters (M, kz, v, E, v, Aoi, N, 
Oexp) take their value directly from the measurements. Three 
adjustable parameters have been systematically varied to choose 
the most relevant value: Cj, has no effect on the results, while 
n and € have been adjusted to reproduce at best the stick- 
slip regime. Ideally, p and a, should not be adjustable, but the 
dispersion in the experimental estimates of œ; is such that their 
values were not sufficiently constrained. In practice, the value 
of p was chosen equal to the one suggested from experiments 
incorporating not only microjunctions within multicontacts, but 
also millimetric smooth sphere/plane individual contacts (Sahli 
et al., 2018). The value of ay was then adjusted to best match the 
overall decay of real contact area during the incipient loading of 
the interface. 

With those values, the time evolution of the tangential load Q 
is quite well-reproduced (see Figure 3). In particular, the slope of 
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FIGURE 4 | Real area of contact vs. tangential force for the six experiments. 
Black (blue) curves show the measured (predicted) data. The red line has slope 
_ and passes through the origin. 
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the incipient loading is correct, which suggests that the stiffness 
used for the individual microjunctions is also correct. In contrast, 
the time evolution of the real contact area is not satisfactory. 
Of course, the total amplitude of the real area decay, from the 
initial contact to macroscopic sliding, is correct, because we start 
from the measured initial value (57;_, A), and we adjusted a, 
to get the correct final value. So we argue that the quality of 
the comparison between the model and experimental results can 
only be assessed through the shape of the real area decay. And as 
can be seen from Figure 4, while the shape of the experimental 
curves A(Q) is essentially quadratic, that of the model curves 
is much more linear [except from the very beginning, when 
all microjunctions are pinned and thus decay quadratically 
according to Equation (3)]. We emphasize that this quasi-linear 
shape is a very robust feature of our model, because we found 
that the predictions are essentially unaffected by changes in the 
model assumptions (elliptic vs. circular microjunctions, Equation 
(3) applied at all times or only before the first depinning event) 
and in the parameter values (for values of 7 and € enabling 
stick-slip or not). 

The shape of the curve A(Q) results from a sum of a 
large number (N) of complex individual behaviors (non-linear 
area decay while pinned, constant area while slipping) with 
distributed parameters (initial area, stiffness, threshold), and 
is therefore unlikely amenable to a simple explanation. We 
can however mention an instructive particular case where 
all microjunctions would have the same initial area. In 
those (unrealistic) conditions, all microjunctions would behave 
identically when submitted to a common displacement X and 
thus depin at the same instant. The total area decay would 
be the sum of N identical quadratic decays, and thus be itself 
quadratic with the total shear load, until macroscopic sliding. 
With those specific (but wrong) initial conditions, we would 


recover a macroscopic area curve with the correct quadratic 
shape and a simple adjustment of the value of a, would allow 
us to provide a good matching with the measurements. This 
example illustrates the major influence of the distribution of 
initial areas on the final shape of A(Q). It also clearly shows 
that the fact that we did not succeed in reproducing a quadratic 
area decay is not a generic problem of our model, but partly 
relates to the initial conditions (through the Ag;) imposed by the 
experimental dataset. 

Could there be other reasons for the failure of our model 
to reproduce the evolution of the real contact area? The main 
model ingredient responsible for this evolution is Equation (3). 
The first possibility is that, despite the evidence brought in Sahli 
et al. (2018, 2019), the anisotropic area reduction under shear 
would not follow a single behavior law at all junction scales, 
from millimeter- to micrometer-sized junctions. This possibility 
is indeed suggested by a recent adhesion-based model of sheared 
sphere/plane junctions (Papangelo et al, 2019), where the 
authors find that the exponent p varies systematically, for a given 
sphere, with the normal load applied, and hence the initial area. 
Here we did not try to apply the model of Papangelo et al. (2019), 
because it would require the knowledge of the characteristic 
radius of curvature of, and normal load on, each individual 
microjunction. In contrast, experimental measurements only 
provide a combination of both quantities, through the area of 
the microjunction. 

We now argue that the solution to the failure of our model will 
presumably be much more complex than a mere improvement 
of the form of Equation (3). The problem may very well be 
that the predicted individual force f; is significantly different 
from the one that really applies on the microjunction. This is 
substantiated by Figure 2 which shows the time evolution of 
the contact area of various microjunctions. Two of them (4 and 
5) were selected to show that the time window over which the 
area decay occurs can be very different: microjunction 4 has 
not started to shrink yet when the decay of microjunction 5 is 
already complete. This observation suggests that the individual 
tangential forces f4 and fs evolve very differently during the 
experiment, although they have very similar initial areas and 
are submitted to the same tangential displacement by the glass 
substrate. We speculate that such a difference may be the result 
of elastic interactions between microjunctions, with junctions in 
a crowed neighborhood evolving differently from those far from 
neighboring junctions!. 

Those interactions are ignored in our model of independent 
microjunctions. We thus believe that, in order for an asperity 
model to have a chance to quantitatively match experiments 
like those considered here, tangential elastic interactions must 
be accounted for to describe the shear behavior of individual 
microjunctions. Such improved models may incorporate those 
tangential interactions in ways similar to models already 
developed for the normal interactions during normal loading of 
rough surfaces (see e.g., Ciavarella et al., 2006; Afferrante et al., 





'The slight initial increase in area of junction 4 in Figure 2 may be due not only to 
elastic interactions but also to a slight aging due to viscous creep. 
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2012) or for the friction of multicontacts (see e.g., Braun and 
Scheibert, 2014; Tromborg et al., 2014; Braun and Peyrard, 2018). 


5. CONCLUSION 


We developed an independent asperity model for the incipient 
shear loading and onset of sliding of dry multicontact interfaces 
between a rough elastic solid and a smooth rigid surface. We used 
it to attempt the first deterministic comparison with experiments 
which, in addition to the macroscopic loads and displacements, 
also considers the individual areas of the many microjunctions 
forming the interface. 

The main outcome is that, although we did our best 
to incorporate experimentally-based behavior laws, parameter 
values and initial conditions into the model, it fails to 
quantitatively reproduce the measurements of Sahli et al. 
(2018, 2019). Based on observations at the microjunction 
scale, we suggest that an interesting starting point for 
future attempts to improve the quantitative deterministic 
comparison between asperity models and experiments, may be 
to incorporate a description of the tangential elastic interactions 
between microjunctions. 

Nevertheless, we anticipate that asperity-based friction 
models, although accounting for tangential elastic interactions, 
may suffer from the same limitations as their counterparts for 
purely normal contact (in particular the difficulty to define 
asperities when a continuum of length scales is involved in 
the topography, see e.g., Miiser et al., 2017; Vakis et al., 2018), 
and may still be unsuccessful to quantitatively match friction 
experiments. We thus urge for the concurrent development of 
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APPENDIX 1: INTEGRATION OF THE 
EQUATION OF MOTION OF THE SLIDER 


Integration Scheme for the Differential 
Equation Giving X(t) 


Equation (1) is integrated by a second order Leapfrog algorithm. 





N 
10 = RO = js | or x) D] nX (A1) 


=l 


is the acceleration of the slider at the beginning of an integration 
step. y(t = 0) has to be calculated at the start of the simulation. 
In our case, we start with X(t = 0) = 0 and X(t = 0) = 0 so 
that y(t = 0) = 0. 

Let dt be the time step. The algorithm first updates X(t) 
according to: 


; i 
X(t + dt) = X(t) + dt X(t) + sae y(t). (A2) 
With this updated value of X we calculate the new acceleration 
y(t + dt). In the second stage of the algorithm X is updated 
according to 


X(t + dt) = X(t) + sd (v(t) + y(t+ dt), (A3) 


and then everything is ready for the next step starting at t + dt. 


To properly select dt we calculate the elasticity constant k; for 
each contact at the start of the simulation. The angular frequency 
of oscillation of the slider of mass M under the total force of the 
contacts is Qo = 4/ Du k;)/M and the corresponding period is 
To = 27 / Qo. The calculation uses dt = To/C+ with C; = 104. 
The exact value of dt depends on the experiment, but a typical 
value is dt = 1 ys. This value is sufficiently small with respect 
to all the frequencies entering in the dynamics of the slider so 
that even a second order algorithm gives a very good numerical 
accuracy. We have however run some calculations with a 4th 
order Runge-Kutta method (Carnahan et al., 1969), which is 
significantly slower, but has errors that decay as dt*, to check the 
accuracy of our calculations. 


Algorithm of the Subroutine Which 


Calculates y 
To compute y(t + dt), X(t + dt) and x; are known. The 
main point is to compute all the forces f; on the junctions. 
The state of each junction is recorded with two flags: 0; 
records its instantaneous state, 0; = 1 for a pinned junction, 
6; = 0 for a slipping junction, and h; keeps track of its 
history, h; = 1 for a junction which has never been slipping 
switches to h; = 0 the first time the junction starts to slip, 
when fi > fi. 

The program scans all the junctions and performs the 
following steps: 


e Compute f; for each junction 


x if hj = 1 the area of the junction depends on fj according 
to Equation (3), which determines €}; and then ki[Ai(fi) 
according to Equation (2). Thus 

fi = kil Af] (X(t + at) — xi) (A4) 
gives an equation for fj. It is too complex to be solved 

analytically. We solve it by an iterative process, using a 

dichotomy method starting from the value of f; from the 

previous step. Once fj is known we update A;(f), k;(A;) and 
fi = 0A; for further steps. 
* if hj =0 
- if 6; = 1 the junction is pinned but A; is fixed, as 
well as k;(A;), and they are known from previous 
iterations so that f; = ki (x (t) — Xi): 
- if0; = 0 the junction is slipping. fi = éfii. 
e Check for transitions in the junction state 
x if 6; = 1 (pinned junction) then if fi > fii the junction 
starts to slip so that 6; switches to 0, f; = ef. hj switches to 0 if 

it was still equal to 1. 

x if 6; = 0 (slipping junction), if IXO < Cmin X 

v the junction repins, 6; switches to 1 and we set x; = 

X — éfa/k;i so that the junction starts in the pinned state 

with fi = efi. 


Once all junctions have been scanned and all f; are determined, 
we can compute y from Equation (A1). 
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